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One of the main diagnostic tools for measuring electron density profiles and the characteristics of long wave- 
length turbulent wave structures in fusion plasmas is Beam Emission Spectroscopy (BES). The increasing 
number of BES systems necessitated an accurate and comprehensive simulation of BES diagnostics, which 
in turn motivated the development of the renate simulation code that is the topic of this paper, renate 
is a modular, fully three-dimensional code incorporating all key features of BES systems from the atomic 
physics to the observation, including an advanced modeling of the optics. Thus renate can be used both 
in the interpretation of measured signals and the development of new BES systems. The most important 
components of the code have been successfully benchmarked against other simulation codes. The primary 
results have been validated against experimental data from the KSTAR tokamak. 
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I. INTRODUCTION 

Although fusion has been an intensely explored area 
for over 50 years, still several key phenomena, such as the 
turbulent transport, spontaneous transition to enhanced 
confinement regimes or the edge localized modes, have 
not been adequately explained. Also a solid predictive ca- 
pability regarding to the confinement properties of fusion 
plasmas is still lacking. The progress in the development 
and testing of theoretical models and predictive simula- 
tion codes strongly relies on measurements from various 
plasma diagnostics. In particular there is a need for ac- 
curate, well localized measurements of density profiles 
and fluctuations, which can be accomplished by Beam 
Emission Spectroscopy (BES). BES is an active plasma 
diagnostic which analyzes the collisionally induced emis- 
sion of either a non-intrusive Diagnostic Neutral Beam 
(DNB) probe or in some cases that of a Neural Beam In- 
jection (NBI) heating bearoi. The emission distribution 
provides information about the distribution of density in 
the plasma. 

Having good spatial resolution (few cms) and high 
sampling rate (few MHz)i BES is ideal for measur- 
ing density fluctuations making it the basic diagnos- 
tic tool for studying the long wavelength fluctuations 
of the plasma turbulence and different transient trans- 
port events. This resulted in a high demand for such 
systems. Since determining the viability of a BES mea- 
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surement conflguration is not trivial, several simulation 
packages have been in use to aid such efforts. The emit- 
ted light intensity can be directly calculated from beam 
and plasma parameters if the cross sections of relevant 
atomic physics processes are provided. A number of 
codes using collisional-radiative models^i^ have been de- 
veloped which can provide the local emission intensity, 
such as FLYCHKi, NOMAD^, SIMULA^, NUBEAMl, 
ALCBEAM^ and the code developed by O. Marchuk^. 
However these arc atomic physics simulation packages 
and accordingly the effects of the observation and beam 
geometry are not taken into account. From these as- 
pects, a more detailed approach is taken in Simula- 
tion of Spectra—, which uses a less precise quasi-steady 
approximation^, but accounts for spectral effects (e.g. 
motional Stark-effect) and includes basic observation and 
device geometry. 

The development of new BES systems and the evalua- 
tion of current measurements required a detailed simula- 
tion which takes all geometrical effects, including obser- 
vation, beam structure and plasma parameter distribu- 
tions into account, allowing it to directly simulate mea- 
sured signals. This motivated the development of the 
renate simulation code. The aim of the current paper 
is to give an overview of renate by providing informa- 
tion about its capabilities and potential applications. 

The remainder of the paper is organized as follows. 
In Section |TT] we describe the main structure of the code 
along with the theoretical models used. Then, in Sec- 
tion lllll an example will be given for the application of re- 
nate along with the validation of renate results against 
a particular experimental setup on KSTAR. Finally, the 
results are discussed in Sec. |TVl 
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II. THE RENATE SIMULATION CODE 

In the present section we give an overview of the BES 
simulation code, renate (Rate Equations for Neutral 
Alkali-beam Technique). The code has a modular struc- 
ture and is written in the IDL languageii. As the name 
suggests, the original purpose of renate was the model- 
ing of alkali atomic beams (lithium and sodium)^^, but 
during the course of its development support for the more 
common hydrogen isotopes was also implemented, re- 
nate has two different atomic physics kernels, both al- 
lowing the modeling of the beam evolution in plasmas 



where Ni is the population of the atomic level i [i = 1 
denotes the ground state and m is the number of con- 
sidered levels), 71/ is the density of species J, R denotes 
the rate coefficients and A denotes the Einstein coeffi- 
cients. The rate coefficients are calculated from the cross 
sections {a) by assuming that the plasma species have a 
Maxwellian population while the beam is assumed to be 
monoenergetic. 

R = {av) = j f{T, v) \v - vi,\ o {\v ~ v^\) dv, (2) 

where vj, is the beam velocity, v is the speed of the par- 
ticle species interacting with the beam, while / (T, v) is 
their Maxwellian velocity distribution. Equation ([1} is 
solved by assuming the plasma parameters to be constant 
during the time needed for the beam to pass the observed 
region. This allows the transformation of Eq. ([T]) into a 
space dependent system of differential equations, which 
can be written in the following matrix form: 

dN 

— ^ \A(x)nJx) + B] • N(x) = C(x) ■ N(x), (3) 
ax 

where the B matrix contains the spontaneous atomic 
transitions while the A(a;) matrix contains all collisional 
processes. A is dependent on the spatial location x 
through the temperature and impurity distribution. Due 
to the aforementioned transformation rule the simulation 
of gas jets is not feasible using Eq. ([3]) as the prescribed 
minimum beam velocity is too high. Eq. ([3]) is solved by 
renate using a fourth order Runge-Kutta scheme. 

In case of hydrogen isotope beams the cross sections for 
(n, I) — >■ (n, Iztl) transitions are so high, that a statistical 
population between the sub-shells can be assumed, this 



with mixed isotope content and arbitrary impurity com- 
position. 

renate is capable of calculating beam evolution by 
using a coUisional-radiative mode]^, taking collisional 
excitation (exc), de-excitation {dexc), ionization (ion), 
charge exchange {CX) reactions and spontaneous de- 
excitation into account. The recombination of the ionized 
beam material and the interaction with the background 
electromagnetic radiation field are neglected. This leads 
to the following time dependent rate equations in the 
frame moving with the beam atoms: 



(1) 
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is called the bundled-n approximation^^ that we adopt 
(here n and I denote the principal and azimuthal quan- 
tum numbers, respectively). Due to the fact that only a 
finite number of atomic levels can be taken into account, 
RENATE solves Eq. (H]) by considering only the first 9 lev- 
els in case of lithium, 8 levels in case of sodium and 6 
levels (shells) in case of hydrogenic species. This approx- 
imation is based on the fact that higher levels tend to be 
only lightly populated (see Fig. [T]) due to the small cross 
sections and high rate of spontaneous de-excitement. Nu- 
merical calculations have shown that increasing the num- 
ber of considered levels above the number considered by 
RENATE only negligibly affects the emission distribution 
(e.g. <1% in case of hydrogen beams). 

The rate coefficients are calculated by the atomic 
physics kernel of renate, using atomic physics data 
from several sources. The cross sections for lithium were 
taken from Schweinzer et alJ^, while data for sodium 
were obtained from Igenbergs et alf^^. Cross section 
data for hydrogenic species were obtained from the IAEA 
ALADDINi^ and the Open ADAS^i databases with the 
corrections from E. Delabie and O. Marchuki^. 

It should be noted that, unlike in the quasi-stationary 
model to be discussed later, the rate coefficients for the 
impurity ions are scaled from the rates of hydrogenic 
plasma species. Both theoretical calculations and exper- 
iments show that the cross section of various processes 
between impurities and plasma species exhibit universal 
behavior which is only dependent on the charge (q) and 
energy (E) of the ioni^ii^. The scaling stipulates that 
the reduced cross section {(Jp'^'^) of process p at the re- 
duced energy [E'^'^'^) is identical for all species. The re- 
duced cross sections and energies for excitation, charge 
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exchange and ionization processes are the following 
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where n^f / = Ry / Ej, is the effective principal quantum 
number with Ry as the Rydberg constant and Ei, as the 
binding energy. 

The rate equation solver of renate has been bench- 
marked against the code developed by O. Marchuk^ for 
hydrogenic beams, and against the SIMULA code^ for 
lithium beams (Fig. [T]), yielding identical results within 
the numerical accuracy of the calculations. 

Instead of solving the rate equations renate can also 
utilize a simpler quasi-stationary model^, taking atomic 
physics data from the Open ADAS databasoii. This 
model, that has low computational requirements, pro- 
vides a quick and approximate solution to the problem 
of the beam evolution. It assumes the life-time of excited 
atomic levels to be negligible compared to the time scale 
of beam evolution. The beam attenuation according to 
the quasi-stationary model is given by 



N{t)=J2N^{t)=No 



r Seff(r)dT 



Seffir) =J2Mr)Sa{na{T),Ta{T),E), (5) 

a 

where N{t) is the total number of non- ionized beam 
particles, Ni{t) is the population of the atomic level i, 
and Sa{na,Ta, E) denotes the beam stopping coefficient 
corresponding to the interaction between the beam and 
plasma species a calculated from rate coefficients^ 



Saina,Ta,E) = 1/ [C" ^ K , T,, i;)] 



(6) 



using the notation of Eq. ([3]) . 

The emission at a given point is taken to be propor- 
tional to the number of excited particles in the relevant 
atomic state and the local density, that is 



A. 
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where A„_i.„' (r) is the emission density corresponding to 
the n — > n' transition at a specific spatial position while 
£„_!.„' (rie, Te, i?) is the effective emission coefficient. 

[Ni'^N^N,) Fi'\ (8) 
( s) 

where An^n' is the Einstein coefficient, Nn is the Saha- 
Boltzmann population for principal quantum shell n and 
f//'' is the effective contribution of ground level excita- 
tions to the population of level n^. 



Both the beam stopping and effective emission coef- 
ficients are provided by the OpenADAS atomic physics 
databaseii. Since the quasi-stationary model neglects 
the finite life-time of excited atomic levels, the results are 
posteriorly smoothed by assuming exponential decay for 
the excited population to account for the effect. In prac- 
tice this means that the emission distribution along the 
beamline is convolved with an exponential decay whose 
characteristic length is set to be the same as the observed 
excited state's. 

The basic assumption of the quasi-static model is that 
the atomic relaxation is much faster than the rate of 
plasma parameter change which is proportional to the 
beam velocity. In some cases this can be violated in the 
plasma pedestal leading to erroneous results (Fig. [5]). 

While there are quite a few simulation codes which can 
calculate beam evolution^i^iiS, some even in SDii^, re- 
nate is also able to take all the subtle geometrical effects 
of the observation into account, renate calculates the 
beam evolution along the beamline in three dimensions, 
while assuming the plasma density, temperature and im- 
purity composition to be fiux functions and distributing 
them according to the provided profiles and magnetic ge- 
ometry of the configuration. 

As it was previously mentioned, EES diagnostic mea- 
surements are performed not only on thin diagnostic 
beams but on heating beams as well. The latter - unlike 
probe beams - are in general not localized to a poloidal 
plane, hence the need for a full 3D simulation. Heating 
beams require large and complex ion sources, thus the 3D 
structure of the neutral beam is not necessarily trivial (as 
shown in Fig. [3]) , which in turn could have a significant 
effect on the measured signal. To account for this, the 
neutral beam itself is modeled as a set of infinitesimally 
thin virtual beams, for which the beam evolution is cal- 
culated individually in the 3D model of the tokamak. 

Having calculated the emission distribution along each 
virtual beam, contributions to the observation channels 
are summed up by the optical modules. There are two 
possible choices for the observation modeling, represent- 
ing different levels of sophistication. The simpler module 
utilizes a pinhole optics approximation which assumes 
that each channel of the optical system can only detect 
light that was emitted within a pyramid whose apex is 
the observation point, see Fig. S) The amount of light 
reaching the detector surface is assumed to be equivalent 
to that reaching the forward aperture, which is calculated 
by integrating the emission density weighted by the geo- 
metrical efficiency factor (coming from aperture distance, 
size and orientation) within the observed volume of space. 

The second, more sophisticated optical module uses 
ray-tracing through the Zema^^^S model of the optical 
system to calculate the observation efficiency for every 
point in space, thus the transfer matrix of the optical 
system (Fig. [5]). This involves a Monte Carlo simulation 
- carried out by Zemax - that creates a large number 
of randomly generated rays originating from a point-like 
source. These are traced through the optical system, tak- 
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FIG. 1: Evolution of the population of atomic levels along the beam normalized to the initial particle number. 
Figure (a) shows the evolution of a lithium beam calculated by RENATE (R) and Simula (S). Figure (b) shows the 
evolution of a hydrogen beam calculated by renate (R) and the code of O. Marchuk (M). 
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FIG. 2: Evolution of the population of the observed 
atomic level along the beam normalized to the initial 
particle number for KSTAR configuration detailed in 
Section Uni 
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FIG. 3: Toroidal projection of the simulated relative 

emission intensity (purple tone color coding) of a 
heating beam with three ion sources (directions with 
blue arrows) in a hypothetical KSTAR configuration, 
also depicting magnetic axis (dotted red) and the edge 
of the plasma (red). 



ing the characteristics of each optical element into ac- 
count. Then the observation efficiency is simply the light 
intensity reaching the individual detector surface divided 
by the emitted intensity. 

This method provides more accurate results, as it takes 
all the optical elements of the system into account, how- 
ever it can only be used in the final design phase of 
BES systems when the detailed optical plan is available. 
Meanwhile the pinhole optics module is much more flexi- 
ble allowing the tryout of a multitude of different config- 
urations without the computationally demanding calcu- 
lation of transfer matrices. It should be noted that while 
RENATE can calculate the Dopplcr-shifted spectrum and 
thus take filtering into account, it is currently unable 
to take other wavelength shifting effects (e.g. motional 



Stark-effect) into account, which means that filtering effi- 
ciency for each channel has to be provided externally (e.g. 
by another simulation code, like Simulation of Spcctrap^). 

The final result of the calculation is the expected ab- 
solute photon currents (1/s) reaching the surface of the 
individual detector segments (Fig. [6]). 

To give a complete account of the capabilities of RE- 
NATE, wc note that it can compute the density pertur- 
bation response matrices, which in turn could be used to 
determine the response in detected signal of the inves- 
tigated BES system to arbitrary density perturbations. 
Using this matrix, the plasma density fluctuations can 
be converted into signal fluctuations, however this is only 
possible if the response is proportional to the strength of 
the perturbation. It is also possible to try and calculate 
the density fluctuations for a given BES system from the 
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FIG. 4: Schematic figure of the pinhole optical model of 
RENATE showing the pyramid shaped observed area of 
space, the beamline and the observation point. 
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FIG. 6: Graphical output of renate showing the 
density along the beamlets (top- red), light profiles 
(top-blue) and the photon currents reaching the 
detector surfaces (bottom) 




FIG. 5: Contour plot of the poloidal projection of the 
transfer matrices of each detector channel created by 
the advanced optical module of renate in a 
hypothetical TEXTOR configuration, also denoting the 
structural elements of the device and the optical paths 
of the terminal line of sights. 



experimental signal data. This aspect of the renate sim- 
ulation code will be thoroughly explored in a follow up 
article. 



III. APPLICATION TO KSTAR TRIAL MEASUREMENT 

The first demonstration of the capabilities of renate 
includes apphcations to the KSTAR BES system^!. With 
the help of RENATE the photon currents reaching the de- 



tector surfaces could be calculated, which - combined 
with the given detector type's characteristics - enabled 
the calculation of the signal-to-noise ratio for each de- 
tector segment, renate was also capable of providing 
estimates of the 3D spatial and 2D poloidal resolution 
of the fluctuation BES system, by calculating the the 
response matrix and integrating it along magnetic field 
lines. The poloidal and toroidal projections of the simu- 
lated configuration are depicted on Fig. [71 

Using measured data the simulations have been re- 
run for the parameters of KSTAR discharge #6123 
at 1744 ms, and its results have been compared 
against experimental results from the KSTAR BES trial 
measurement^^. The magnetic geometry was recon- 
structed using EFIT— while the temperature profile was 
obtained from Electron Cyclotron Emission (ECE) mea- 
surements, thus assuming T"; = Tp (Fig. [8]). Since at the 
time of the measurements KSTAR had only a single in- 
terferometry channel for plasma density measurement, it 
was not possible to get localized density data from the ex- 
periment. Instead a profile of specific shape was fitted to 
the line integrated data. This density profile is fiat in the 
core plasma, while follows the shape of the relative BES 
emission in the edge as shown in Fig. [5] Impurity data 
were also unavailable thus a homogeneous 1% carbon im- 
purity was assumed, corresponding to Zeff ~ 1.35. 

During the experiment only the central beam source 
was active, resulting in a single deuterium beam (13 A, 
90 keV) with known current distribution. The measure- 
ment was carried out using a 4 x 8 matrix of detectors 
with spatial resolution of 1 cm. Subtracting the back- 
ground light from the measured data and taking into 
account the transmittance of the optical system along 
with the effect of the filtering and detector amplification 
yielded the absolute photon current reaching the aper- 
ture for each individual channels. This method however, 
carries a significant relative error mainly because only a 
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FIG. 7: Toroidal and poloidal projections of the simulated configuration for KSTAR BES, showing some structural 
elements (black), the flux surfaces (blue), the magnetic axis (red circle), the relative emission density (purple tone 
color coding) and the edge lines of sights (red) from the observation point 
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FIG. 8: Plasma parameter profiles used in the BES 
modeling of the KSTAR BES trial measurement. 
Density data (blue) were fitted to interferometry results 
while temperature data (red) were obtained from ECE 
measurements. 



relative calibration of the detector amplification was pos- 
sible from gas shots. By estimating this error term to be 
roughly 30% and assuming the other terms to be indepen- 
dent, the measured data have a total 33% relative error. 
These data can be compared against RENATE results as 
shown in Fig. (S) The density profile of this simulation 
was highly arbitrary and a basic pinhole optics module 



FIG. 9: Measured and simulated photon currents for 
individual detector channels for KSTAR #6123 at 1744 
ms. The detectors are sorted into a 4 x 8 matrix, 
indexing start from the detector of the first row 
observing the innermost part of the plasma. The error 
estimates are high due to inaccurate input data and 
experimental uncertainties. 



was used to model the observation, which represent error 
terms of 40% and 20% respectively. This means an esti- 
mated relative error of 44% for the the simulation results. 

As shown in Fig. IH] the simulated signal is close to 
the measurement. The upper part of Fig. [TU] shows the 
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FIG. 10: Normalized light profiles of the individual 
detector rows (upper figure) and the relative error of the 
simulation compared to the measurement (lower figure). 



normalized light profiles of the individual detector rows, 
where it is evident that the shape of the simulated signal 
matches that of the measurement. This is no surprise as 
the density profile was fitted to the relative BES emis- 
sion which is roughly proportional to the local density as 
shown by Equation ([3]). The real test of the code is the 
lower part of Fig. [TUl which shows the relative error of 
the simulation compared to the measurement. 

On average the simulated results arc 60% higher than 
the measurement. Considering that the standard devia- 
tion of the error is rather small (20%) the dominant term 
is caused by a systematic error which can be caused by 
a number of effects. 

• The density profile is the main source of error as it 
is the primary input of RENATE simulations. Since 
only line integrated density was available, the shape 
of the density profile in the non-BES region of the 
plasma could cause a systematic error. Using a 
parabolic profile for the plasma center instead of 
the flat one shown in Fig. [8] could mitigate this ef- 
fect. It was found that a 40% reduction of the BES 
region density eliminates this term leaving only an 
error of only 20%, which could easily be attributed 
to other uncertainties (e.g. pinhole optics model). 

• The RENATE simulation utilized a simple pinhole 
optical model which could have easily neglected a 
portion of the real observed area. 

• The atomic physics data of most hydrogenic tran- 
sitions have a relative error of 10-20 %^^i^^ . 



IV. CONCLUSIONS 

This paper presents an overview of the renate code, 
which is designed to help the development of new BES 
systems and the evaluation of experimental data. It is 
a modular BES simulation capable of simulating most 
key levels of a BES system. It possesses two different 
types of atomic physics kernel, allowing it to calculate 
beam evolution either by solving the rate equations or by 
employing a quasi-steady approximation. Both of these 
methods utilize the latest atomic physics data, renate 
is able to model the three-dimensional geometry of the 
simulated system in detail, including the inner structure 
of the beam and the optical setup, allowing it to model 
complete BES systems, in particular heating beam mea- 
surements. 

The rate equation solver module of renate has been 
benchmarked for hydrogenic species and lithium. Mean- 
while the experimental results for the KSTAR BES trial 
measurement agree with the simulated results well within 
the estimated error. In the future renate is to be com- 
pared against other experimental results and is to be up- 
graded to be able to take more wavelength shifting effects 
(e.g. motional Stark-effect) into account. 
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